{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Constraint Optimization\n",
    "\n",
    "## Objective and Prerequisites\n",
    "\n",
    "If you are looking to improve your modeling skills, then try this tricky constraint optimization problem. We’ll show you how to model this problem as a linear programming problem using the Gurobi Python API and solve it using the Gurobi Optimizer.\n",
    "\n",
    "This model is example 18 from the fifth edition of Model Building in Mathematical Programming by H. Paul Williams on pages 273 and 328-330.\n",
    "\n",
    "This modeling example is at the advanced level, where we assume that you know Python and the Gurobi Python API and that you have advanced knowledge of building mathematical optimization models. Typically, the objective function and/or constraints of these examples are complex or require advanced features of the Gurobi Python API.\n",
    "\n",
    "**Download the Repository** <br /> \n",
    "You can download the repository containing this and other examples by clicking [here](https://github.com/Gurobi/modeling-examples/archive/master.zip). \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Problem Description\n",
    "\n",
    "In an integer programming model, the following constraint occurs:\n",
    "\n",
    "$$\n",
    "9x_{1} + 13x_{2} -14x_{3} + 17x_{4} + 13x_{5} - 19x_{6} + 23x_{7} + 21x_{8} \\leq 37\n",
    "$$\n",
    "\n",
    "All the decision variables in the constraint are binary, and the objective is to find another constraint involving the same binary variables that is logically equivalent to the original constraint, but that has the smallest possible absolute value of the right-hand side (RHS).\n",
    "\n",
    "It is well established in the theory of  integer programming that the ability to efficiently solve an  integer programming problem depends on how \"close\" the linear programming relaxation is to the convex hull of all the integer solutions of the \n",
    "integer programming problem. A linear programming relaxation is \"tight\" when this formulation is \"close\" to the convex hull of all the integer solutions. We say that an equivalent constraint with respect to an original constraint is \"tighter\" whenever the equivalent constraint has the same integer solutions as the original one, but has removed some fractional feasible solutions of the original one. These ideas are illustrated in the following figure.\n",
    "\n",
    "![OptimizedConstraint](OptimizedConstraint.PNG)\n",
    "\n",
    "The original constraint in the figure is determined by the line AB and the \"tighter\" equivalent constraint is defined by the line CD. Notice that the equivalent constraint has the same set of integer solutions and has removed fractional feasible solutions of the original constraint.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Model Formulation\n",
    "\n",
    "A procedure for simplifying a single 0–1 constraint has been described by Bradley et al. (1974). We adopt their procedure of using a linear programming model. It is convenient to consider the constraint in a standard form with positive coefficients in descending order of magnitude. This can be achieved by the transformation:\n",
    "\n",
    "$$\n",
    "y_{1} = x_{7}, y_{2} = x_{8}, y_{3} = 1 - x_{6}, y_{4} = x_{4},\n",
    "$$\n",
    "\n",
    "$$\n",
    "y_{5} = 1 - x_{3}, y_{6} = x_{5}, y_{7} = x_{2}, y_{8} = x_{1}.\n",
    "$$\n",
    "\n",
    "The resulting equivalent constraint is\n",
    "\n",
    "$$\n",
    "23y_{1} + 21y_{2} + 19y_{3} + 17y_{4} + 14y_{5} + 13y_{6} + 13y_{7} + 9y_{8} \\leq 70.\n",
    "$$\n",
    "\n",
    "All inequalities of the form $\\sum_{i=1}^{8} b_{i}*y_{i} \\leq b_{0} $ equivalent to the constraint to be optimized are characterized by a set of inequalities of a linear programming model defined in terms of the decision variables $b_{i}$. Observe that the variables $b_{i}$ are the coefficients of the transformed constraint.\n",
    "\n",
    "In order to capture the total logical import of the original constraint, we search for subsets of the indices known as *roofs* and *ceilings*. Ceilings are ‘maximal’ subsets of the indices of the variables for which the sum of the corresponding coefficients does not exceed the right-hand side coefficient. Such a subset is maximal in the sense that no subset properly containing it, or to the left in the implied lexicographical ordering, can also be a ceiling. For example, the subset {1, 2, 4, 8} is a ceiling, $23 +21 + 17 + 9 \\leq 70$, but any subset properly containing it (e.g. {1, 2, 4, 7, 8}) or to the ‘left’ of it (e.g. {1, 2, 4, 7}) is not a ceiling. Roofs are ‘minimal’ subsets of the indices for which the sum of the corresponding coefficients exceeds the right-hand side coefficient. Such a subset is ‘minimal’ in the same sense as a subset is ‘maximal’. For example {2, 3, 4, 5} is a roof, \n",
    "$21 +19+ 17 + 14 > 70 $ ,\n",
    "but any subset properly contained in it (e.g. {3, 4, 5}) or to the ‘right’ of it (e.g.\n",
    "{2, 3, 4, 6}) is not a roof. See the details about how these constraints are derived in Bradley et al. (1974).\n",
    "\n",
    "Therefore, all inequalities of the form $\\sum_{i=1}^{8} b_{i}*y_{i} \\leq b_{0} $ equivalent to the constraint to be optimized are characterized by the following constraints.\n",
    "\n",
    "**Ceiling constraints**: <br />\n",
    "\n",
    "$$\n",
    "\\{1,2,3 \\}: b_{1} + b_{2} + b_{3} \\leq b_{0}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\{1,2,4,8 \\}: b_{1} + b_{2} + b_{4} + b_{8}  \\leq b_{0}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\{1,2,6,7 \\}: b_{1} + b_{2} + b_{6} + b_{7}  \\leq b_{0}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\{1,3,5,6 \\}: b_{1} + b_{3} + b_{5} + b_{6}  \\leq b_{0}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\{2,3,4,6 \\}: b_{2} + b_{3} + b_{4} + b_{6}  \\leq b_{0}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\{2,5,6,7,8 \\}: b_{2} + b_{5} + b_{6} + b_{7} + b_{8} \\leq b_{0}\n",
    "$$\n",
    "\n",
    "**Roof constraints**: <br />\n",
    "\n",
    "$$\n",
    "\\{1,2,3,8 \\}: b_{1} + b_{2} + b_{3} + b_{8} \\geq b_{0} + 1\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\{1,2,5,7 \\}: b_{1} + b_{2} + b_{5} + b_{7}  \\geq b_{0} + 1\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\{1,3,4,7 \\}: b_{1} + b_{3} + b_{4} + b_{7} \\geq b_{0} + 1\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\{1,5,6,7,8 \\}: b_{1} + b_{5} + b_{6} + b_{7} + b_{8} \\geq b_{0} + 1\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\{2,3,4,5 \\}: b_{2} + b_{3} + b_{4} + b_{5} \\geq b_{0} + 1\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\{3,4,6,7,8 \\}: b_{3} + b_{4} + b_{6} + b_{7} + b_{8} \\geq b_{0} + 1\n",
    "$$\n",
    "\n",
    "\n",
    "$$\n",
    "\\text{Decreasing value:} \\quad b_{i} \\geq b_{i+1} \\quad \\forall i=1..7\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\text{Non-negativity:} \\quad b_{i} \\geq 0 \\quad \\forall i=0,1..7\n",
    "$$\n",
    "\n",
    "The objective is to minimize the RHS of the original constraint.\n",
    "$$\n",
    "\\text{Objective function:} \\quad \\text{Min} \\;\\; (b_{0} - b_{3} - b_{5} )\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Python Implementation\n",
    "\n",
    "We import the Gurobi Python Module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install gurobipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gurobipy as gp\n",
    "from gurobipy import GRB\n",
    "\n",
    "# tested with Python 3.11 & Gurobi 11.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# list of indices of variables b[i]\n",
    "\n",
    "indices = [*range(9)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Deployment\n",
    "We create a model and the decision variables. These variables determine the coefficients  of the transformed equivalent constraint that minimized the value of the RHS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using license file c:\\gurobi\\gurobi.lic\n"
     ]
    }
   ],
   "source": [
    "model = gp.Model('ConstraintOptimization')\n",
    "\n",
    "# Variables b[i]\n",
    "b = model.addVars(indices, name=\"b\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ceiling constraints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Constraint {1, 2, 3}:\n",
    "\n",
    "C123 = model.addConstr( b[1] + b[2] + b[3] <= b[0] , name='C123')\n",
    "\n",
    "# Constraint {1, 2, 4, 8}}:\n",
    "\n",
    "C1248 = model.addConstr( b[1] + b[2] + b[4] + b[8] <= b[0] , name='C1248')\n",
    "\n",
    "# Constraint {1, 2, 6, 7}:\n",
    "\n",
    "C1267 = model.addConstr( b[1] + b[2] + b[6] + b[7] <= b[0] , name='C1267')\n",
    "\n",
    "# Constraint {1, 3, 5, 6}:\n",
    "\n",
    "C1356 = model.addConstr( b[1] + b[3] + b[5] + b[6] <= b[0] , name='C1356')\n",
    "\n",
    "# Constraint {2, 3, 4, 6}:\n",
    "\n",
    "C2346 = model.addConstr( b[2] + b[3] + b[4] + b[6] <= b[0] , name='C2346')\n",
    "\n",
    "# Constraint {2, 5, 6, 7, 8}:\n",
    "\n",
    "C25678 = model.addConstr( b[2] + b[5] + b[6] + b[7] + b[8] <= b[0] , name='C25678')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Roof constraints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Constraint {1, 2, 3, 8}:\n",
    "\n",
    "R1238 = model.addConstr( b[1] + b[2] + b[3] + b[8] >= b[0] + 1 , name='R1238')\n",
    "\n",
    "# Constraint {1, 2, 5, 7}:\n",
    "\n",
    "R1257 = model.addConstr( b[1] + b[2] + b[5] + b[7] >= b[0] + 1 , name='R1257')\n",
    "\n",
    "# Constraint {1, 3, 4, 7}:\n",
    "\n",
    "R1347 = model.addConstr( b[1] + b[3] + b[4] + b[7] >= b[0] + 1 , name='R1347')\n",
    "\n",
    "# Constraint {1, 5, 6, 7, 8}:\n",
    "\n",
    "R15678 = model.addConstr( b[1] + b[5] + b[6] + b[7] + b[8] >= b[0] + 1 , name='R15678')\n",
    "\n",
    "# Constraint {2, 3, 4, 5}:\n",
    "\n",
    "R2345 = model.addConstr( b[2] + b[3] + b[4] + b[5] >= b[0] + 1 , name='R2345')\n",
    "\n",
    "# Constraint {3, 4, 6, 7, 8}:\n",
    "\n",
    "R34678 = model.addConstr( b[3] + b[4] + b[6] + b[7] + b[8] >= b[0] + 1 , name='R34678')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Decreasing value coefficients\n",
    "These constraints guarantee the ordering of the coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Decreasing value of coefficients\n",
    "\n",
    "DV1_2 = model.addConstr( b[1] >= b[2] , name='DV1_2')\n",
    "DV2_3 = model.addConstr( b[2] >= b[3] , name='DV2_3')\n",
    "DV3_4 = model.addConstr( b[3] >= b[4] , name='DV3_4')\n",
    "DV4_5 = model.addConstr( b[4] >= b[5] , name='DV4_5')\n",
    "DV5_6 = model.addConstr( b[5] >= b[6] , name='DV5_6')\n",
    "DV6_7 = model.addConstr( b[6] >= b[7] , name='DV6_7')\n",
    "DV7_8 = model.addConstr( b[7] >= b[8] , name='DV7_8')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The objective is to minimize the  value of the RHS side of the original constraint."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Objective function\n",
    "\n",
    "model.setObjective(b[0] -b[3] -b[5] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gurobi Optimizer version 9.1.0 build v9.1.0rc0 (win64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 19 rows, 9 columns and 76 nonzeros\n",
      "Model fingerprint: 0x88e2e477\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+00]\n",
      "  Objective range  [1e+00, 1e+00]\n",
      "  Bounds range     [0e+00, 0e+00]\n",
      "  RHS range        [1e+00, 1e+00]\n",
      "Presolve time: 0.00s\n",
      "Presolved: 19 rows, 9 columns, 76 nonzeros\n",
      "\n",
      "Iteration    Objective       Primal Inf.    Dual Inf.      Time\n",
      "       0   -2.0000000e+30   7.000000e+30   2.000000e+00      0s\n",
      "      16    2.5000000e+01   0.000000e+00   0.000000e+00      0s\n",
      "\n",
      "Solved in 16 iterations and 0.01 seconds\n",
      "Optimal objective  2.500000000e+01\n"
     ]
    }
   ],
   "source": [
    "# Verify model formulation\n",
    "\n",
    "model.write('OptimizeConstraint.lp')\n",
    "\n",
    "# Run optimization engine\n",
    "\n",
    "model.optimize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "________________________________________\n",
      "The equivalent optimized constraint is:\n",
      "\n",
      "6.0x[1] + 9.0x[2] + -10.0x[3] + 12.0x[4] + 9.0x[5] + -13.0x[6]  + 16.0x[7] + 14.0x[8] <= 25.0\n"
     ]
    }
   ],
   "source": [
    "# Output report\n",
    "\n",
    "\n",
    "# Optimized values of the coefficients resepect to the original variables x.\n",
    "\n",
    "a = {}\n",
    "\n",
    "a[0] = b[0].x - b[3].x - b[5].x\n",
    "a[1] = b[8].x\n",
    "a[2] = b[7].x\n",
    "a[3] = -b[5].x\n",
    "a[4] = b[4].x\n",
    "a[5] = b[6].x\n",
    "a[6] = -b[3].x\n",
    "a[7] = b[1].x\n",
    "a[8] = b[2].x\n",
    "\n",
    "print(\"________________________________________\")\n",
    "print(f\"The equivalent optimized constraint is:\")\n",
    "\n",
    "\n",
    "optimized_constraint = {}\n",
    "for i in indices:\n",
    "    if i==0:\n",
    "        optimized_constraint[i]=f\"{a[i]}\"\n",
    "    else:\n",
    "        optimized_constraint[i]= f\"{a[i]}x[{i}]\"\n",
    "        #print(optimized_constraint)\n",
    "            \n",
    "print(f\"\\n{optimized_constraint[1]} + {optimized_constraint[2]} + {optimized_constraint[3]} + {optimized_constraint[4]} + {optimized_constraint[5]} + {optimized_constraint[6]}  + {optimized_constraint[7]} + {optimized_constraint[8]} <= {optimized_constraint[0]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## References\n",
    "\n",
    "H. Paul Williams, Model Building in Mathematical Programming, fifth edition.\n",
    "\n",
    "Bradley, G.H., Hammer, P.L. and Wolsey, L. (1974) Coefficient Reduction for Inequalities in 0–1 Variables. Mathematical Programming, 7, 263–282.\n",
    "\n",
    "Copyright © 2020 Gurobi Optimization, LLC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
